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Full quantum statistical NVT simulation of the five-particle system Hj^ has been carried out 
using the path integral Monte Carlo method. Structure and energetics is evaluated as a function 
of temperature up to the thermal dissociation limit. The weakly density dependent dissociation 
temperature is found to be around 4000 K. Contributions from the quantum dynamics and thermal 
motion are sorted out by comparing differences between simulations with quantum and classical 
nuclei. The essential role of the quantum description of the protons is established. 
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The triatomic molecular ion H;^ 



is a five-body system 
consisting of three protons and two electrons. Being the 
simplest polyatomic molecule it has been the subject of 
a number of theoretical and experimental studies over 
the yearsir— . Experimentally, the H;^ ion was first de- 
tected in 1911 by ThompsonS, however, definite spectro- 
scopic studies were carried out not until 1980 by Okapi. 
Since then, this five-body system has proven to be rele- 
vant, also in astrophysical studies concerning the inter- 
stellar media and the atmosphere of gas planets. There- 
fore, low-density high-temperature Hg ion containing at- 
mospheres have been studied experimentally^ as well as 
computationally^. 

Until now, the computational approaches have consis- 
tently aimed at finding ever more accurate potential en- 
ergy surfaces (PES) for \i^ at zero Kelvin, and conse- 
quent calculations of the rovibrational states^"'^^. These 
calculations include Born-Oppenheimer (BO) electronic 
energies in various geometries often supplemented with 
adiabatic and relativistic correctionsiSJ^. For the study 
of rovibrational transitions it is desirable to have an ana- 
lytical expression for the PES, which is usually generated 
using Morse polynomial fitsiS. Inclusion of the nonadia- 
batic effects, however, has turned out to be a cumbersome 
task, and so far, they have not been rigorously taken into 
account^.. 

In this work, we evaluate the full five-body quantum 
statistics of the Hg ion in a stationary state at temper- 
atures below the thermal dissociation at about 4000 K. 
We use the path integral Monte Carlo (PIMC) approach, 
which allows us to include the Coulomb correlations be- 
tween the particles exactly in a transparent way. Thus, 
we are able to monitor the fully nonadiabatic correlated 
quantum distributions of particles and related energies 
as a function of temperature. Furthermore, we are able 
to model the nuclei as classical mass points, in thermal 
motion or fixed as conventionally in quantum chemistry, 
and find the difference between these and the quantum 
delocalized nuclei. 

The PIMC method is computationally expensive, but 
within the chosen models and numerical approximations 
it has been proven to be useful with exact correlations 
and finite temperature^"—. For zero Kelvin data with 



benchmark accuracies, however, the conventional quan- 
tum chemistry or other Monte Carlo methods, such as the 
diffusion Monte Carlo^^, are more appropriate. Thus, it 
should be emphasized that we do not aim at competing 
in precision or number of decimals with the other ap- 
proaches. Instead, we will concentrate on physical phe- 
nomena behind the finite-temperature quantum statis- 
tics. 

Next, we will briefly describe the basics of the PIMC 
method and the model we use for the ion. In the results 
and discussion section we first compare our 160 K PIMC 
"ground state" to the zero Kelvin ground state, and then, 
consider the higher temperature effects. 



II. METHOD 

According to the Feynman formulation of the quantum 
statistical mechanics^ the partition function for interact- 
ing distinguishable particles is given by the trace of the 
density matrix: 
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S is the action, /? — l/k-oT, 
T — /3/M, Rm = Rq and M is called the Trotter num- 
ber. In this paper, we use the pair approximation in the 
action— iSi for the Coulomb interaction of charges. Sam- 
pling in the configuration space is carried out using the 
Metropolis procedure^^ with bisection movea^^. The to- 
tal energy is calculated using the virial estimator^^. 

The error estimate in the PIMC scheme is commonly 
given in powers of the imaginary time time-step t^ 
Therefore, in order to systematically determine thermal 
effects on the system we have carried out all the simula- 
tions with T = 0.03£'jJ^, where E^i denotes the unit of 
Hartree. Thus, the temperatures and Trotter number M 
become fixed by the relation T — (fceMr)"^. 

In the following we mainly use the atomic units, where 
the lengths, energies and masses are given in units of the 
Bohr radius (ao), Hartree (-Eh) and free electron mass 
{me), respectively. 

The statistical standard error of the mean (SEM) with 
2SEM limits is used as an error estimate for the observ- 
ables, unless otherwise mentioned. 



III. MODELS 

Two of the five particles composing the Hg ion are 
electrons. For these, we do not need to sample the exact 
Fermion statistics, but it is sufficient to assign spin-up 
to one electron and spin-down to the other one. This is 
accurate enough, as long as the thermal energy is well 
below that of the lowest electronic triplet excitation. 

We do the same approximation for the three protons, 
too. This is even more safe, because the overlap of well lo- 
calized nuclear wave functions is negligible and related ef- 
fects become very hard to evaluate, anyway. On the other 
hand, however, the nuclear exchange due to the molecu- 
lar rotation results in the so called zero-point rotations. 
These too contribute to energetics less than the statis- 
tical accuracy of our simulations. Therefore, we ignore 
the difference between ortho-H;^ (/ = 3/2) and para-H;^ 
(J = 1/2). Thus, the protons are modeled as "boltz- 
mannons" with the mass nip = 1.83615267248 x lO'^mg. 
The higher the temperature, the better is the Boltzmann 
statistics in describing the ensemble composed of ortho- 
and para-H|^. 

For the NVT simulations we place one H^ ion into a 
cubic box with the volume of (SOOao)'^ and apply pe- 
riodic boundary conditions (PBC) and minimum im- 
age principle. This corresponds to the mass density of 
^ 1.255 X 10~^ gcm~'^. This has no essential effect at 
low T, but at high T the finite density gives rise to the 
molecular recombination balancing the possible dissocia- 
tion. Within the considered temperature range the dis- 
sociations are very rare. 

The electrons are always simulated with the full quan- 
tum dynamics. For the nuclei, however, we use three 
models to trace the quantum and thermal fluctuations, 
separately. The case of full quantum dynamics of all par- 
ticles we denote by AQ (all-quantum), the mass point 
model of protons by CN (classical nuclei) and the adi- 
abatic case of fixed nuclei by BO (Born-Oppenheimer 
potential energy surface). 



IV. RESULTS AND DISCUSSION 
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Figure 1; (Color online) Total energy of the H^ molecular 
ion as a function of temperature. Fully nonadiabatic quan- 
tum statistical simulations, AQ (blue circles), classical nuclei 
simulations, CN (red squares), and the equilibrium geometry 
Born-Oppenheimer simulation, BO (black triangles). Zero 
Kelvin datai^i^^ii is given for comparison: BO ground state en- 
ergy at equilibrium internuclear geometry (black dash-dotted 
line), energy including the nuclear zero-point motion (green 
dashed line) and energy at the barrier to linearity (grey solid 
line). 2SEM statistical error estimate is shown by the er- 
ror bars from simulations at the H;^ 
~ 1.255 X 10"'' gcm-^ 



ion density (300oo) or 



energy. Note however, that the nuclear spins and zero 
point rotation are not included in our model of Hg . 

The lowest electronic excitation from the BO ground 
state is a direct Franck-Condon one (O.TlOi^H)*''^ to dis- 
sociative potential curve: H;^ ^-112 + H^ or HJ -^ 
H^ + H.--'— The dissociation energies (Dg) are 0.169ii^H 
and 0.241£'h, respectively. 

The linear geometry with equal bond lengths 1.53912ao 
{Dooh) is a saddle point on the BO PES at 
-1.27868190£^Hii or 0.06515£;h above the BO energy at 
the equilibrium geometry. This energy is usually called 
as the barrier to linearity^. The zero Kelvin energetics is 
shown in Fig. [1] by the three horizontal lines. 



A. Ground state: zero Kelvin reference data 

The equilibrium geometry of the Hg ion in its ground 
state is an equilateral triangle D^h for which the inter- 
nuclear equilibrium distance is i? = 1.65ar/. The best 
upper bound for the electronic ground state BO energy 
to date is -1.34383562502£:k^. The vibrational normal 
modes of H;^ are the symmetric-stretch mode i^i and 
the doubly degenerate bending mode 1^2 ■ The latter one 
breaks the full symmetry of the molecule, and therefore, 
it is infrared active^. 

The vibrational zero-point energy is 0.01987i?H, and 
the so called rotational zero-point energies are 0.00029£'h 
and 0.00040i?H for para- and ortho-H;^, respectively^^. 
These yield about 0.020215i?H for the average zero-point 



B. PIMC ground state: 160 K 

At our lowest simulation temperature, T « 160 K, 
the electronic system is essentially in its ground state. 
For the total energy we find - 1.3438 (2) Sh, see the BO 
black triangles in Fig. [T] The thermal energy is ksT = 
0.000507i?H, and therefore, the contribution from the ro- 
tational and vibrational excited states is also small and 
we find — 1.3406(29)i?H, see the CN red square in same 
Fig. The full quantum simulation includes vibrational 
zero-point contribution and yields — 1.3233(12)£'h, about 
0.0205(14)iJH above the BO energy in a good agreement 
with about 0.0202£;h in Refs.^'". 

From our AQ simulation we still find the equilateral 
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Figure 2: (Color online) Nuclear pair correlation functions 
(bond length distributions) at different temperatures from the 
quantum statistical simulations (solid lines), and from the 
classical nuclei simulations (dashed lines). The zero Kelvin 
equilibrium internuclear distance is given as a vertical black 
dash-dotted line. The distributions include the r^ weight and 
normalization to unity. (Note that the r^ weight is usually 
not included in description of extended or periodic systems) 
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Figure 3: (Color online) Expectation values of the internu- 
clear distance at different temperatures from distributions in 
Fig. [2] Quantum statistical simulations (blue circles) and clas- 
sical nuclei simulations (red squares). The FWHM limits are 
shown by triangles (all the lines are for guiding the eye). The 
zero Kelvin equilibrium internuclear distance is shown as a 
horizontal black dash-dotted line. 



triangle configuration of the nuclei with the internuclear 
distances increased to (R) = 1.723(4)ao, which indicates 
an increase of about 0.073(4)ao, as compared with the 
zero Kelvin BO equilibrium distance bond lengths. In- 
terestingly, within the error limits this is the same as the 
bond length increase of the hydrogen molecule ion HJ. 
The zero-point energy of Hj^ is about 2.7 times as large 
as that of the H^ ion—, as expected from the increase 




r (units of a ) 



Figure 4: (Color online) Proton-electron pair correlation 
functions at the four temperatures from the full quantum sta- 
tistical simulations, AQ (solid lines), and from simulations 
with the classical nuclei, CN (dashed lines). That from the 
BO scheme is given at the lowest (electronic) temperature, 
only (dash-dotted line). Notations are the same as in Fig. [J] 
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Figure 5: (Color online) Electron-electron pair correlation 
functions from the same simulations as those in Fig. U Nota- 
tions are the same as in Figs. [2] and (4] 



of vibrational modes from one to three — the zero-point 
energy of our model does not contain the rotational zero- 
point energy, as mentioned earlier. 

The thermal motion (CN), alone, increases the bond 
length to (R) = 1.658(4)ao, only, see the data in Figs. [5] 
and Fig. [3] This clearly points out the difference between 
quantum and thermal delocalization of nuclei at low T. 

For the proton-electron and electron-electron inter- 
actions the differences between our two approaches are 
smaller than in the proton-proton case but still distinc- 
tive. Comparison of the fixed nuclei simulation to the 
CN one shows that the two schemes give almost identical 
distributions. The AQ distributions, however, cannot be 
labeled identical with those from the CN or fixed nuclei 



simulations. The distributions are given in Figs. |3] and [SI 
where the notations are the same as in Fig. [51 

The calculations of the relativistic corrections involve, 
among other things, evaluation of the contact densities, 
{6{rij)), for the electron-nuclei and the electron-electron 
pairsi^. For the electron-nuclei contact density at the 
BO equilibrium configuration we get 0.1814(20), and for 
the AQ case 0.1765(20). For the electron-electron pair 
we get 0.0182(3) and 0.0166(3), for BO and AQ ap- 
proaches respectively. The estimated uncertainties due 
to extrapolation to the contact are given in parenthe- 
sis. The zero Kelvin reference valueai^ for the BO case 
are 0.181242 (electron-nuclei) and 0.01838663 (electron- 
electron). Thus, the quantum dynamics of the nuclei 
turns out to be significant factor in lowering the contact 
densities, too. 

See the snapshot of the AQ simulation in Fig. [5| for 
some intuition of the low-temperature quantum distribu- 
tions in imaginary time. 



C. High temperature phenomena 
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Figure 6: (Color online) xy-plane (z-projection) snapshot of 
the Hg" ion from quantum statistical simulation with Trotter 
number 2^^, i.e. temperature of about 160 K, for all particles. 
"Polymer rings" describing the electrons are in the background 
(green and blue) and those of the nuclei are placed on top 
(yellow) . 



With the increasing temperature the increasing con- 
tribution from rovibrational excitations is clearly seen in 
the total energies shown in Fig.[T] Contributions from the 
electronic excitations do not appear, because the lowest 
excitation energy 0.710£'h is much too high as compared 
to the thermal energy fc^T. Consequently, the equilib- 
rium geometry BO energy depends on the temperature 
almost negligibly. For convenience, the essential energet- 
ics related data has been collected into Table [H also. 

As expected, the increase in the total energy due to the 
classical rovibrational degrees of freedom is 9 x ^fc^T, 
defining the slope of the CN line. The most prominent 
quantum feature in AQ curve is, of course, the zero-point 
vibration energy. At higher temperatures, however, by 
comparing the AQ and CN curves we see that the quan- 



Table I: Energetics of the HJ molecular ion. The energies 
are given in the units of Hartree (atomic units). Simulation 
data is given with 2SEM error estimates. BO refers to Born- 
Oppenheimer calculation at equilibrium geometry. The refer- 
ence data is rounded to convenient accuracy. The "barrier to 
linearity" is 0.06515-Bh ^ 1.8 eV above the Ebo at K. 
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tum nature of nuclear dynamics becomes less important, 
except for dissociation. 

At the dissociation limit we find the molecule with 
quantum nuclei somewhat more stable than the one with 
classical nuclei. With the relatively low density, (300ao)^, 
the molecule is mainly kept in one piece above 4000 K in 
the former case, whereas more dissociated in the latter. 
The total energy becomes higher for the CN than the 
AQ case slightly below 4000 K, see Table [D The total 
energies at this crossing point are above the "barrier to 
linearity" )Siii already. 

At higher temperatures, T > 4100 K, other configura- 
tions, such as H2 -I- H"*", HJ -|- H and 2H -I- H"*", start play- 
ing more significant role in the equilibrium dissociation- 
recombination processes. These will considered in our 
next study. 

The nuclear pair correlation function or bond length 
distributions. Figs. [2] and [3j follow the energetics dis- 
cussed, above. There, the zero-point vibration in AQ 
case is seen even better. At the zero Kelvin limit both 
the expectation value and the distribution, in particular, 
are significantly different from those of the CN case. 

The temperature dependence in the other pair correla- 
tion functions is weak, see Figs.[3|and[Sl Obviously, this is 
the case, because electrons do not present a quantum-to- 
classical transition in the temperature range considered, 
now. Thus, the evolution in distributions in Figs.[3|and[5| 
following the rising temperature arises from the changes 
in the nuclear dynamics, and mostly, from the change in 
the conformation or the bond lengths, presented in Fig. [31 



V. CONCLUSIONS 

In this study, the path integral Monte Carlo method 
was shown to be a successful approach for examination of 
quantum statistics of the five-particle molecule, H^ ion. 
The method is based on the finite temperature mixed 
state description, and thus, it gives information, which 
is complementary to the high-accuracy zero Kelvin de- 
scription of conventional quantum chemistry. It was also 
shown how contributions from quantum and thermal dy- 
namics to particle distributions and correlation functions 
can be sorted out, and furthermore, quantum to classical 
dynamics transition can be monitored. 

Our approach is fully basis set and trial wave function 
free. It is based on the Coulomb interactions, only, and 
allows the most transparent interpretation of consequent 
particle-particle correlations. 

Simulation at 160 K essentially reproduces the zero 
Kelvin data from conventional quantum chemistry. Of 
course, a proper extrapolation to K can be done for 
more accuracy. Born-Oppenheimer (BO) potential en- 
ergy surface and the equilibrium geometry can be found 
by using classical nuclei with fixed coordinates. De- 
scription of the zero-point motion within our nonadia- 
batic five-body quantum simulation gives the vibrational 
zero-point energy accurately. We find an increase of 
0.073(4)ao in the bond length due to the nonadiabatic 
zero-point vibration. The classical thermal contribution 
at 160 K is 0.008(4)ao, only. 



With the raising temperature the rovibrational excita- 
tions contribute to the energetics, as expected, whereas 
the electronic part remains in its ground state in the spirit 
of BO approximation. At about 4000 K the H^ ion disso- 
ciates, weakly depending on the ion density. We find that 
the full quantum molecule dissociates at slightly higher 
temperature compared to the one, where the nuclei are 
modeled by classical particles with thermal dynamics, 
only. Thus, we conclude the necessity of the quantum 
character of the protons in the correct description of dis- 
sociation. 

We find that the nuclear quantum dynamics has a dis- 
tinctive effect on the pair correlation functions, too. This 
is least for the electron-electron pair correlation func- 
tion, stronger for the electron-proton one and largely 
increased in the proton-proton correlations. These are 
seen in the contact densities, and consequently, in the 
relativistic corrections where relevant. 
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